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Non-Euclidean plates are thin elastic bodies having no stress-free configuration, hence exhibiting 
residual stresses in the absence of external constraints. These bodies are endowed with a three- 
dimensional reference metric, which may not necessarily be immersible in physical space. Here, based 
on a recently developed theory for such bodies, we characterize the transition from flat to buckled 
equilibrium configurations at a critical value of the plate thickness. Depending of the reference 
metric, the buckling transition may be either continuous or discontinuous. In the infinitely thin 
plate limit, under the assumption that a limiting configuration exists, we show that the limit is a 
configuration that minimizes the bending content, amongst all configurations with zero stretching 
content (isometric immersions of the mid-surface). For small but finite plate thickness we show the 
formation of a boundary layer, whose size scales with the square root of the plate thickness, and whose 
shape is determined by a balance between stretching and bending energies. 

PACS numbers: 46.25.Cc, 46.70.De 87.10.Pq 



I. INTRODUCTION 

The classical literature on thin elastic bodies deals pri- 
marily with two types of bodies — plates and shells. Math- 
ematically, a plate can be viewed as a continuous stack 
of identical flat surfaces glued together, whereas a shell 
can be viewed as a continuous stack of non-identical (and 
not necessarily flat) surfaces glued together. The term 
non-Euclidean plate was coined in [1] to describe thin 
elastic bodies which, like plates, do not exhibit struc- 
tural variations across their thin dimension, and yet, 
unlike plates, do not have a planar rest configuration. 
Such elastic bodies can neither be described as shells, 
which bear structural variations across their thin dimen- 
sion (e.g., shells do not display reflectional symmetry 
about the mid-surface), and possess curved stress-free 
rest configurations. Non-Euclidean plates exhibit resid- 
ual stresses even in the absence of external constraints, 
and are therefore inherently frustrated. 

Elastic bodies having such properties are ubiquitous 
in biology. Growing tissues, such as plant leaves, are 
relatively thin elastic structures that may exhibit com- 
plex equilibrium configurations in the absence of exter- 
nal forces [2, 3]. In addition, thin elastic bodies that have 
no stress-free configuration have been engineered in a 
laboratory [4], for example, the environmentally respon- 
sive gels shown in Figure 1 . 

There are various ways to treat elastic bodies which 
exhibit residual stress. One way is to treat the residual 
stress as a physical field and characterize its properties 
[5]. Another is to decompose the deformation gradient 
into a product of a plastic (or growth) process, which de- 
forms the body from a rest configuration into some "vir- 
tual" configuration, and an elastic relaxation from the 
virtual configuration to the current configuration [6, 7] . 
A third approach is to decompose the strain tensor addi- 




FIG. 1: (Color online) Four elastic plates made of thermo- 
responsive gel as described in [4]. All four structures bear 
no structural variation across their thickness. Their radius is 
3 cm. The mid-surface of the positively curved discs (a and 
b) possess a reference metric of constant Gaussian curvature 
K = 0.11 cm"^ . The mid-surface of the negatively curved sur- 
faces (c and d) possess a reference metric of constant Gaussian 
curvature of opposite sign K = -0.11 cm'^. Plates (a) and (c) 
are 0.75 mm thick, whereas plates (b) and (d) are 0.6 mm thick. 



tively into a plastic (or growth) strain, leading to a virtual 
configuration, and an elastic strain [8, 9]. The treatment 
presented here for residually stressed three-dimensional 
bodies is very similar to the third approach. All defor- 
mations which are not of elastic nature are completely 
ignored, i.e. it is assumed that the virtual configura- 
tion which the growth or plastic deformation led to is 
known, and the appropriate elastic relaxation is solved. 
This in turn enables us to treat large "plastic strains" in 
a non-iterative manner. 

A static theory of non-Euclidean plates was developed 
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in [1 ] following the fundamental principles laid by Trues- 
dell [10] and its modem interpretation by Ciarlet [11, 12]. 
The starting point in [1] is the formulation of a covari- 
ant three-dimensional elasticity theory in the form of 
an energy functional. A first notable property of this 
energy functional is its expression in terms of the three- 
dimensional metric of the configiiration. Specifically, the 
energy density is quadratic in the deviation of the met- 
ric from a reference metric. This deviation of the metric 
is a strain, which reduces to the standard Green-Saint 
Venant strain for bodies that have a rest configuration. 
The second notable property of our model is that the 
reference metric is not required to have a vanishing Rie- 
mannian curvature tensor, i.e., it may not be immersible 
in (hence the name "incompatible elasticity theory"). 
As a result, there exist no rest configurations in which 
the strains vanish ever5rwhere, hence the state of frus- 
tration. In a second step, a two-dimensional elasticity 
theory is derived, using a generalization of the standard 
Kirchhoff-Love assumptions [13, 14]. The end result is 
an energy functional which depends on surface proper- 
ties of the midplane of the plate, namely, on the first and 
second fundamental forms. Like in the classical Foppl- 
von-Karman theory [15] the energy functional is a sum of 
a stretching term and a bending term. The bending term 
is minimized in flat configurations, whereas the stretch- 
ing term measures, in an sense, deviations of the 2D 
surface metric from a prescribed reference metric, and 
vanishes only in surfaces which are isometric immer- 
sions of the given 2D metric. The lack of iiranersibility 
of the three-dimensional metric manifests in the lack of 
a planar stretching-free configuration. 

At this stage we have a new model, which we be- 
lieve to be applicable to a large variety of physical and 
biological systems, whose properties are governed by 
essentially two-dimensional shaping mechanisms. In 
[1], a single application was demonstrated for the case 
of an unconstrained thin plate, whose two-dimensional 
reference metric is that of a punctured spherical cap. A 
buckling transition was shown to occur at a critical plate 
thickness. 

In this paper we study two behaviors exhibited by 
unconstrained non-Euclidean plates. First, we study 
the transitions from flat to buckled equilibrium states 
as the plate thickness crosses a critical value — ^the buck- 
ling threshold. We derive an explicit expression for the 
critical thickness in terms of the stress field in the pla- 
nar configuration. An immediate implication is that the 
plane-stress solution always becomes unstable for suf- 
ficiently thin plates, provided that it is not trivial, i.e., 
that the stress is not identically zero. We apply this 
analysis to three reference geometries of constant Gaus- 
sian curvature of different types — an elliptic, a flat and 
a hyperbolic metric. We show that the buckling tran- 
sition may be either continuous (super-critical) or dis- 
continuous (sub-critical). In particular, we show that 
the buckling threshold may deviate significantly from 
the so-called crossover point, which is based on a bal- 



ance between the plane-stress energy and the energy that 
minimizes the Willmore functional. We show an exam- 
ple in which the crossover thickness underestimates the 
buckling threshold by more than an order of magnitude. 

Second, we analyze the equilibrium configurations 
and energies in the limit where the plate thickness tends 
to zero. We show that if a limit configiiration exists, 
then it is the minimizer of the bending content amongst 
all configurations with zero stretching content, i.e., the 
Willmore energy minimizer amongst all isometric im- 
mersion of the 2D reference metric [16]. For a small 
but finite thickness, deviations from isometric immer- 
sions are more pronounced near the free boundary of 
the domain, forming a boundary layer, which we obtain 
in explicit form. In particular, the size of this boundary 
layer is foimd to scale with the square root of the plate 
thickness. 



II. THEORY OF NON-EUCLIDEAN PLATES 

In this section we briefly review the modeling of non- 
Euclidean plates, first described in [1 ] . The starting point 
is a three-dimensional covariant elasticity theory, based 
on the principles of hyper-elasticity [10]: the elastic en- 
ergy is a volume integral over an energy density, which 
depends only on (i) the local value of the metric ten- 
sor and (ii) local characteristics of the material that are 
independent of the configuration (the use of the metric 
tensor, rather than the deformation, as primitive vari- 
able has been originally proposed by Antman [17], and 
has been recently advocated by Ciarlet and co-workers 
[11, 18-21]). 

Let Q c be an elastic body endowed with a set 
of material curvilinear coordinates x = {x^,x^,x^) c 
D c R^. Let r denote the mapping from the do- 
main of parametrization D into Q — r{x) is called the 
configuration — then the induced Euclidean metric is 
gij = djr ■ djT, where di = dldx\ Our model assumes 
the existence of a reference metric gij{x), such that the 
elastic energy density vanishes at a point x if and only if 
the actual metric coincides with the reference metric at 
that point, gij{x) = gij{x). While the reference metric is 
required to satisfy the properties of a metric — it is S5an- 
metric positive-definite — it is not necessarily immersible 
in R'', hence the name of the theory as "incompatible 
three-dimensional elasticity" 

The strain tensor is defined as half the deviation of the 
metric from the reference metric, 

^ij = \iSii - gij)- 

It coincides with the Green-Saint Venant strain tensor in 

the case where there exists a rest configuration, and the 
curvilinear coordinates form a Cartesian parametriza- 
tion in the rest configuration, i.e., when gij = 5ij. For 
small deviations of the metric from the reference metric, 
the energy functional is truncated at the first non-trivial 
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term, i.e., it is quadratic in the strain tensor, yielding 

E= r w{g)^/\g\dx^dx^dx^, (1) 



where 



and 



wig) = -A''"eije,i, 



A'V* = Ar¥' + f^(rV + f'¥')' (2) 
where A, are the Lame coefficients. 
Comments: 

1. We adopt the Einstein summation convention 
whereby repeated indices imply summation. 

2. Latin lowercase characters i, • • • = 1, 2, 3 are used 
to denote indices of three-dimensional tensors. We 
will use below Greek characters a,|S, ■■■ - 1,2 to 
denote indices of two-dimensional tensors. For 
any tensor fl;y, \a\ denotes its determinant. 

3. The tensor g'> is the tensor reciprocal to gij. The 
raising and lowering of indices is only defined with 
respect to the reference metric. For example, the 
tensor g^' is defined as g^^g'^gu and not as the re- 
ciprocal of gij, which we denote by 

4. The volume element in (1) is determined by the ref- 
erence metric rather than the actual metric. Note 
that it is not a priori clear whether the volume ele- 
ment should be derived from the reference metric 
or from the actual metric. In any case, the differ- 
ence between the two choices is of higher order in 
the strain. 

5. The structure (2) of the elastic tensor is imposed by 
the assumption of spatial isotropy. 

6. In standard (or "compatible") nonlinear elasticity, 
the energy density is sometimes written in terms 
of the Euclidean distance of the deformation gradi- 
ent, Vr, from the group of proper rotations, SO(3). 
In the same spirit, the energy density in an incom- 
patible elasticity theory may be expressed in terms 
of 

dist(Vr, .^g), 

where is the set of matrices, R, such that R^R = 

7. In summary, given the reference metric g, the elas- 
tic problem is formulated as follows: find the met- 
ric g that minimizes the energy functional (1), sub- 
ject to the constraint that it is embeddable, in par- 
ticular, that the corresponding Riemann curvature 
tensor vanishes. 



With a three-dimensional elasticity theory in hand, we 
focus the attention on plate-like structures. We define 
a plate to be an elastic body for which there exists a 
parametrization in which the reference metric takes the 
form 



gn gi2 O"* 

g21 g22 

1 



(3) 



with gij independent of x^. The plate is called even if D = 
S X [-t/2, t/2] with § c and t a constant, and thin if t 
is much smaller than all other dimensions. We identify 
the two-dimensional tensor ga^ as the metric tensor of 
a siirface. By assirmption it is constant across the plate 
thickness. It is easy to see that the three-dimensional 
reference metric (3) is immersible in R'' if and only if the 
two-dimensional reference metric gap has zero Gaussian 
curvature. 

To derive a reduced two-dimensional energy density 
in terms of the mid-surface configuration we used an 
adaptation of the Kirchhoff-Love assumptions [13, 14]: 
We first assume that the stress is parallel to the mid- 
surface, and then that £,3 = (the order in which 
the two assumptions are used is essential). Integrat- 
ing the energy functional (1) over the thin direction, it 
takes straightforward manipulations to derive an en- 
ergy functional which depends on the first and sec- 
ond fundamental forms of the mid-surface. Specifically, 
let f{x^,x^) = r{x^,x^,0) be the immersion of the mid- 
surface, then 

gaj} = daf ■ dfif and ha^ = dj^daf ■ 3\f, 

are the first and second fundamental forms, where N is 
the unit vector normal to the mid-surface. The energy 
functional is given by 



W = tEs + t^EB, 



(4) 



where 



Es= I wsdS Eb= \ 
Js Js 



wb dS 



are called the stretching and bending contents, 

WS = lA"l'y'{gap-ga^){gy6-gy6) 

are their respective densities, where 



(5) 



j^afiyb 



V 



2(1 +v) \1 - V 



2v 



and dS = ^J\g\dx^dx^ is the infinitesimal surface element. 
The coefficients Y and v are Yoimg's modulus and the 
Poisson ratio, which can be related to the Lame coeffi- 
cients. The elastic energy is positive definite for constant 
values of Y > and -1 < v < 1. 

Comments: 
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1. Like in the classical Foppl-von-Karman and Koi- 
ter theories, the energy functional is a sum of (i) 
a stretching energy, which scales with the plate 
thickness and attains a minimum for an isomet- 
ric immersion of the mid-plane surface, and (ii) a 
bending energy, which scales like the third power 
of the thickness and attains a minimum for flat 
configurations. The equilibriirai configuration is 
the minimizer of the sum of both stretching and 
bending energies. 

2. Rather than working with the energy W given 
by (4) we will work with the energy-per-unit- 
thickness, 

W o 

E = — = Es + fEs. (6) 

Henceforth, we will refer to E as "the en- 
ergy". Thus the stretching energy is i-independent 
whereas the bending energy scales with P. Obvi- 
ously, both W and E have the same minimizer. In 
addition, we rescale the units of energy by a factor 
of Y/(l -I- v), such that the tensor A"^y^ takes the 
final form 

3. A different derivation of an elastic functional sim- 
ilar to (6) may be foimd in [12]. Eq. (6) may be 
identified as the elastic energy in the Koiter shell 
model when the "target bending tensor" ha^ is set 
to zero [22]. 

4. With the above rescaling the stretching and bend- 
ing density contents can be written in the more 
compact form 

= 8(i3^[tr(g-v - + 1 m-'g - D'] 

5. The energy functional is expressed in terms of 
the first two fundamental forms ga^i and haji of 
the mid-surface. The two forms are not indepen- 
dent: they must satisfy the three Gauss-Mainardi- 
Codazzi compatibility conditions, 

where the Christoffel symbols are given by 



6. The two-dimensional stress and moment tensors 
are defined as 

s"^ = A^^y^eys and m"^ = ^A"^y%6, 
so that 

1 1 

Ws + fwB = ^S^l^eaj} + ^m"<^hai}. 

7. A surface f{x^,x^) will be called an isometric im- 
mersion if the two-dimensional metric gap coin- 
cides with the two-dimensional reference metric 
ga^, i.e., if the stretching energy is zero. In the 
case of an isometric immersion the bending con- 
tent density wb can be identified with the density 
of the Willmore functional, 

^^=24(1^^^)' 

where H and K are the mean and Gaussian curva- 
tures of the surface [16]. Note that since X is an 
isometric invariant, its value is prescribed by the 
reference metric. 

8. In summary, the two-dimensional elastic problem 
is defined as follows: given the two-dimensional 
reference metric gap, find a symmetric positive 
definite tensor field gap, and a symmetric tensor 
field hap, that together minimize the energy func- 
tional (6), subject to the constraint that the Gauss- 
Mainardi-Codazzi equations are satisfied. 

III. THE INFINITELY-THIN PLATE LIMIT 

In many applications, the elastic body is thin to an 
extent that the equilibrium configuration (of its mid- 
surface) remains practically tmchanged upon further 
thinning. In other words, we identify an asymptotic 
regime, which we may call the infinitely-thin plate limit, 
which in oui model corresponds to the limit f — > 0. It 
is important to stress that there are also opposite cases, 
where the thinner the body is, the more convoluted the 
equilibrium configuration is, with no evidence that a 
t — > limit exists (e.g., in [23] a torn plastic sheet exhibits 
a self -similar shape, whose cut-off scale is comparable to 
the thickness of the sheet). 

Under the assumption that g admits an isometric 
embedding of finite bending content, and that a t — > 
(weak) limit configuration exists, we show in Ap- 
pendix A that the limit is a minimizer of the Willmore 
fimctional amongst all isometric embeddings. 

The first assumption that the bending content is finite 
is non-trivial. If it does not hold, then a limit configu- 
ration may not exist. We expect, however, the second 
assumption, regarding the existence of a (weak) limit, to 
become eventually superfluous, yet further analysis is 
required before this assumption can be relaxed. 
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IV. THE BUCKLING TRANSITION 

A. Plane-stress solution 

A configuration is a flat surface if ha^ = 0. The con- 
figuration f{x^, x^) that minimizes the elastic energy un- 
der the constraint that the surface be flat is called the 
plane-stress solution. It is the minimizer of the stretching 
content, which is given by 

E = l ^/^^a^dS = \ j^^'^^'^'igafi - ga^)igy6 ' gy6)dS, 

with respect to all flat metrics ga^. To find the minimizer, 
we consider an in-plane perturbation of the stirface. 

The reason we perturb the configuration rather than the 
metric is that the three components of the configtiration 
are independent, whereas the three entries of the metric 

tensor are constrained by the Gauss-Mainardi-Codazzi 
relations. The corresponding variation of the metric is 

Using the fact that d^daf ■ dyf = T^^ggriy, the energy vari- 
ation is 



hence A"Py^ = l{g"^gl^^ + g^^gf^'')- Denote = r and 
x^ = 6 and consider a reference metric in semi-geodesic 
parametrization of the form. 



with 0(r) yet to be specified. The domain of parametriza- 
tion is 

(r,0)e[rmin,w]x [0,271), 

with periodicity in the 0-axis, so that the topology of the 
body is that of a punctured disc. 

The equilibrium configuration is expected to preserve 
the axisymmetry of the intrinsic geometry, hence we seek 
plane-stress solutions of the form 

f{r, 9) = {(p{r) cos 9,(p{r) sin 0,0). 

Elementary calculations show that the resulting two- 
dimensional metric is 

_/(c/)')2 0\ 

where cp' = dcp/dr, from which we derive the two- 
dimensional stress tensor. 



5E = ^ s'^^[g„y{dpvy) + Tlg,„vy]dS + 0{v'). ' igy6-gy6)g - 2[ q (cp^/^ -D/O^j- 

(8) 



Integrating by parts, requiring the first variation to 
vanish for any perturbation v'^ (and using the identity 
dagfiy = r^jsgriy " ^lygr,?)' the Euler-Lagrange equations 



Finally, the Christoffel symbols are given by 



are 



r - r/<t>' 

-(p/(p' 



and r - 1 ° '^'/'^ 

(9) 



(7) hence the resulting plane-stress equation is 



with boundary conditions s"^n^, = 0, where Up is the 
outward unit vector tangent to the plate and normal to 
its boimdary. We refer to (7) as the plane-stress membrane 
equations. Note that the plane-stress equations do not 
depend on the plate thickness t, which only comes into 
play when there is a competition between stretching and 
bending energies. 

Comment Equation (7) is the Euler-Lagrange equa- 
tion associated with the energy functional (6) when 
t = Q. It is expected to hold in the limit t — > even 
for non-flat configurations. In general, (7) constitutes 
two equations for three unknown functions (the three 
components of the metric tensor gap), i.e., the system is 
under-determined. In the present case, the third equa- 
tion which removes this under-determinacy is that the 
curvatures be identically zero. 

Examples In the following examples we consider for 
simplicity the case of a vanishing poisson ratio (v = 0), 



;(cl)<//[(f)2-l]) 



1 , 



(10) 



with boundary conditions (p'{rjnm) = (p'i^max) = 1- 

We solve the plane-stress equation (10) for three fami- 
lies of metrics: 



1. A family of elliptic metrics, 

<^(r) = —— sin ^^Kr, 



(11) 



where K > is the constant Gaussian curvature 
of the reference metric. Although such a metric is 
consistent with an infinite set of immersions, the 
immersion that minimizes the Willmore functional 
is easily identified — it is a (punctured) spherical 
cap. 



2. A family of conical flat metrics, 

0(r) = ar, 



(12) 
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with a < 1. Here the isometric immersion that 
minimizes the Willmore functional has the form of 
a truncated cone (a circular frustum). Note that 
although the reference metric is flat, all isomet- 
ric immersions have non-zero bending energy due 
to the topological constraint (periodicity in the 6 
axis). 

3. A family of hyperbolic metrics, 

0(r) = ^— sinh V^r, (13) 

where X < is the constant Gaussian curvature. 
Unlike the two former cases, the minimizer of the 
Willmore fimctional amongst all isometric embed- 
dings is not known explicitly, yet, it is known that 
isometric embeddings with finite bending content 
do exist [24]. 

The plane-stress solutions are shown in Figures 2-4 for 
the domain 0.1 <r< 1.1. The solutions were obtained by 
a simple shooting procedure, with a fourth-order adap- 
tive ODE solver. For each metric we plot the solution 
(p{r) along with the spatial profile of the stress compo- 
nents sl{r) and Sg(r) given by (8) up to the lowering of 
one index (the reason for displaying stress components 
with mixed upper and lower index is that only then all 
the components have the same units, hence can be com- 
pared). In the three cases, the solution (p{r) is close to 
linear. Significant differences are however observed in 
the stress components. 

For the elliptic metric (Figure 2) the body is in a state 
of compression along the r direction (s[; < 0), whereas 
it is compressed in the 6 direction near the inner radius 
and stretched along the 6 direction near the outer ra- 
dius. For the flat metric (Figiire 3) the situation is similar 
with compression ever5where along the radial direction, 
while the angular stress switches from extension in the 
vicinity of the inner boundary to compression at larger 
radii, and again extension in the vicinity of the outer 
boundary. Finally, for the hyperbolic metric (Figure 4) 
the radial stress is ever5rwhere positive (i.e., in a state 
of extension), whereas the angular stress is in a state of 
extension very close to the inner radius and in a state 
of compression at large enough distances form the cen- 
ter. In Figure 5 we show two toy models generating 
hyperbolic and elliptic geometries, which elucidate the 
behavior of the azimuthal (hoop) stress. 

B. Stability analysis and buckling threshold 

Let f{x^,x^) be the plane-stress configuration. Any 
small enough perturbation can be decomposed into a 
sum of in-plane and out-of-plane displacements, 

where "N is the unit vector normal to the surface and dyf 
is the covariant derivative defined in (B2). Given such 




r 



FIG. 2: (a) Plane-stress solution cj)(r) for the elliptic metric (11) 
with Gaussian curvature K = 1. (b)-(c) The corresponding 
principal stresses s].{r) and Sg{r). 




r 



FIG. 3: (a) Plane-stress solution (p(r) for the flat metric (12) with 
a = 0.58. (b)-(c) The corresponding principal stresses s'^r) and 

a perturbation we calculate in Appendix B the variation 
in the elastic energy, 

6E = J" (dws + 6wb^ dS, 

where the variation in stretching content density, 6ws, 
is given by (B6) and the variation in bending content 
density, Swg, is given by (B16); for flat surfaces these ex- 
pression simplify considerably as hap = 0. Note that the 
plane-stress solution enters in the energy variation both 
through the stress s"^ and through the metric parameters, 
gap and Tl^. 

The plane-stress solution is locally stable if the energy 
variation is positive for every choice of sufficiently small 
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FIG. 4: (a) Plane-stress solution (p{r) for the hyperbolic metric 
(13) with Gaussian curvature K = -1. (b)-(c) The correspond- 
ing principal stresses sj:(r) and Sg{r). 




(b) 




FIG. 5: Cartoons of hyperbolic (a) and elliptic (b) plates. In 
both cases, two punctured discs are "glued" one inside the 
other and forced to remain planar. In the hyperbolic case, the 
inner perimeter of the outer disc is too long, compared with the 
outer perimeter of the inner disc. In such a case the irmer disc 
is stretched azimuthally, while the outer disc is compressed 
azimuthally. In the elliptic case, the inner perimeter of the 
outer disc is too short, hence the inner disc is compressed az- 
imuthally, while the outer disc is stretched azimuthally. The 
color bar on the right represents the azimuthal strain at equi- 
librium (computed numerically). A similar cartoon was first 
presented in [25]. 



(non-trivial) perturbation. As is well-known, local sta- 
bility can be determined by considering only the leading 
order terms (in powers of v, w) in the energy variation. 
The defining property of the plane-stress solution is that 
the terms that are linear in the in-plane perturbation c^' 

(the integral of 60;^^'"^ in (B6)) vanish for every choice of 
c^'. Thus, to leading order, the energy variation decom- 
poses into a sum of terms that are quadratic in v and 
terms that are quadratic in w, 



6Ef°\v) 



6Ef^\w)= f ls"l^{Vaw){Vpw)dS 
Js 2 

J^_Lyi«ft'<5(v^^V„zi;)(V6Vya;)rfS 



(14) 



(the subscripts {i, j) refer to the power of v and w). 

Since, to leading order, the energy variation is de- 
composed into a sum of a c-dependent term and a zv- 
dependent term, the minimization can be performed on 
each component separately. By assumption, the plane- 
stress solution is the energy minimizer with respect to 
in-plane perturbations, thus minimum energy variation 
is obtained for v'!' =0. 

It remains to consider the energy variation due to out- 
of-plane perturbations. The bending term 6Eg''^\w) is 
always positive due to the positive-definiteness of the 
tensor yL''^'* [39]. Whether the stretching term 6Ef'^\w) 
is sign-definite depends on the plane-stress solution. In 
fact, if the stress tensor is not everywhere positive-definite, 
then there exists a perturbation w for which bE^^'^\io) is 
negative, and by taking the plate thickness f sufficiently 
small, the total energy variation can be made negative. 
We have thus recovered the following general result: 

Given a reference metric, the plane-stress solu- 
tion is linearly stable against buckling, indepen- 
dently of the plate thickness, only if the stress is 
everywhere positive-definite. In other words, an 
infinitely thin plate cannot sustain compression 
without buckling. 

We will now show that the existence of a buckling 
threshold is always guaranteed, unless the plane-stress 
solution is trivial, i.e., s"'^ = (which in turn occurs only 
if the reference metric is flat). We start by noting that the 
plane-stress equations (7) can be rewritten as 



= 0, 



where V^- is the covariant derivative defined in the ap- 
pendix. Let now x{^^>^^) be a scalar field satisfying 
V^VaX = and consider the integral 



J^s«/'(V„x)(V^X)'^S. 



bE = 6Ef°\v) ■ 



The surface element dS is defined in terms of the ref- 
erence metric. Writing dS = ( y/lgl/ y/lgl) yjlgldx^dx^, we 
may now integrate by parts (the covariant derivative 



[6Ef^\w) + t^ 6Ef'^\iv)] + 0{v^,v^w,vw'^,w%a.tislies the usual rules of integration by parts provided 
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that the surface element is consistent with the Christoffel 
symbols), using the boundary conditions s"l^n^ = 0, 



X 



Since the covariant derivative satisfies the Leibniz rule 
for the derivative of products, it follows from the plane- 
stress equations and the definition of x that / = 0. 

Thus, if there exists a scalar function x that has a non- 
zero (covariant) gradient and satisfies ^a^^X = 0/ then 
the fact that / = implies that s"^ is not ever5rwhere 
positive-definite. A simple way to show that such a func- 
tion does exist is to endow the planar equilibrium state 
with a Cartesian set of coordinates. Then the covariant 
derivative reduces into a simple partial derivative, and 
the function, say, xi^^rX^) = has the desired property. 

We may summarize as follows: 

A sufficiently thin unconstrained non-Euclidean 
plate will always buckle unless the plane-stress 
solution is trivial, i.e., s^f' = 0. 

Equation (14) provides a characterization of the critical 
thickness t = th at which buckling first occurs. At criti- 
cality, t = tj,, there exists a non-trivial (i.e., non-uniform) 
perturbation which to leading order is marginally unsta- 
ble, i.e., 

inf r {ls^^V„w){V^w)+!^A"^y\V^VaW){V6VyW)]dS 
which implies that 



By the above analysis this supremum is guaranteed to 
be non-negative, and zero if and only if s"^ =0. 

Comments: 

1. Equation (15) provides a mean for generating 
lower bounds for the buckling threshold by choos- 
ing appropriate trial functions, w. 

2. The energy variation (14) is a quadratic functional 
of w of the form. 



6Ef^\w) + t^5Ef^\w) = {w,:Kw), 



(16) 



= 0, 



where (•, •) is the standard irmer-product on S and 
is a self-adjoint second-order differential oper- 
ator. Above the buckling threshold 'K is positive- 
definite. The buckling threshold ti, corresponds to 
the largest t for which has a zero eigenvalue. 
From a numerical point of view, the latter char- 
acterization is the easier way for computing the 
buckling threshold. 
Examples We turn back to the punctured discs con- 
sidered in the previous subsection. Note that in all three 
cases there exist negative stress components, hence a 
buckling transition is guaranteed to occur at some finite 
thickness. 

We denote by (p{r) the solution to the plane-stress equa- 
tion (10). For an out-of-plane perturbation w{r, 6), sub- 
stituting (9), we get 



-12 J^s"l\Vaiv){Vpw)dS 
wl7^t^t /g.A«^>'*(V^V„K;)(VaVyU;)dS" 



i= sup 



(15) 



(V^V^i.) = idadpw)-Tl^id,w) = ^^^(^/ J). ^ lct>{w'/4>') J 

where we denote by prunes derivatives with respect to 
r and by dots derivatives with respect to 6. Thus, 



^ Jo Jr„ 



s"{w'f + s^''^Udrdd 



bEf\w) 



mm 

•In 



"24 Jo i„ 



w 



+ 1 



(D2 



(17) 



with s"" and s^^ given by (8). Due to the periodicity in 6 
it is natural to expand the perturbation in Fourier series, 

OO CO 

w{r,d) = ao{r) + V2^a„(r) cosn0 + V2^fc„(r) sinn0. 

n=l n=l 

Because (17) is quadratic in w, both terms reduce into a 
sum over Fourier components, 

OO 

5Ef^\w) = 2^[6E^(a„) + dE"gib„)] 

n=0 

OO 

6Ef^\w) = J][6E^(fl„) + 6E'^{b„)], 

n=0 



where we define for every function z = z(r). 
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6E^(z) = - {s"'{z'f + s^\nzf] O dr 

-| Z^'' max I 



r 



nz 



(54 



0.2 0.4 0.6 0.8 1.0 1.2 1.4 



0.164 0.233 0.285 0.329 0.367 0.401 0.431 



TABLE I: Buckling threshold i\, versus the Gaussian curvature 
K for the elliptic geometry (11). 



The buckling threshold (15) is given by 



tl = sup 



aM Ln=omiia„) + 6El{b„)] 



Corollary 1 Let 



(Q^ = sup 



-6e;;(z) 



Then it is clearly the case that < ti, for every n. On the 
other hand, tj, < max„ , which together implies that 

*^="^rT^£«(^- 

Thus, tmless the buckling transition is degenerate, then 
the marginally stable perturbation at the bifurcation 
point involves a single Fourier mode. 

Corollary 2 A buckling transition occurs if either 

or s"" are somewhere negative. Suppose that s''''(r) > 0, 
i.e., the radial stress is ever5rwhere extensional. It follows 
that 6E°(z) > for every z (every axisymmetric pertur- 
bation, n = 0, increases the stretching energy). If s'^" is 
somewhere negative, then there exist non-axisymmetric 
perturbations that reduce the elastic energy. That is, the 
buckling transition breaks the axial symmetry. 

Elliptic geometry Consider first the elliptic geometry 
(11) for the same parameters as in Figure 2. The buckling 
threshold occurs at t{, = 0.367, and corresponds to an 
axisymmetric mode (n = 0). The critical mode is shown 
in Figure 6a. In Table I we show the buckling threshold 
tt, versus the Gaussian curvature K. As expected, the 
buckling threshold is higher the more curved the surface 
is. 

Flat geometry Consider next the flat geometry (12) 
for the same parameters as in Figure 3. The buckling 



o 



(a) 



o 

T3 
O 

E 



O 



(b) 




-1 



FIG. 6: Critical modes for the elliptic (a), flat (b) and hyperbolic 

(c) geometries. For both elliptic and flat geometries the critical 
mode is axisymmetric (n = 0), hence we only display a cross 
section. In the hyperboKc case the first mode to destabilize is 
n = 3. 



-K 



0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 



0.0768 0.110 0.135 0.157 0.175 0.192 0.208 0.222 0.236 0.248 



TABLE II: Buckling threshold h versus the Gaussian curvature 
-K for the hyperbolic geometry (13). In all cases the critical 
mode has harmonic n = 3. 



threshold occurs at = 0.387, also for an axis5rmmetric 
mode. The critical mode is shown in Figure 6b. 



Hyperbolic geometry Consider finally the h5rper- 
bolic geometry (13) for the same parameters as in Fig- 
ure 4. Since s" > it follows that the critical mode 
must break the polar symmetry. Indeed, the least stable 
mode, which changes stability at th = 0.1845, has har- 
monic n = 3. It is depicted in Figure 6c. Note how lower 
is the buckling threshold for the hyperbolic geometry. 
Finally, we show in Table II the buckling threshold 
versus the Gaussian curvature K. 
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C. Buckling threshold versus crossover point 

Equation (15) expresses the buckling threshold ti, as 
a supremum over trial normal deflections. As such, it 
provides an easy way to generate lower bounds for the 
buckling threshold. An approximation often used to es- 
timate the buckling threshold is the so-called crossover 
point between the lowest-energy isometric immersion 
and the plane-stress solution. In this section we show 
that the crossover point can often yield a significant un- 
derestimate to the buckling transition. 

The equilibrium configuration is the one that mini- 
mizes the energy ftmctional (6). Two upper botmds that 
correspond to extreme cases are the plane-stress solution, 
which involves zero bending energy, and the isomet- 
ric immersion that minimizes the Willmore functional, 
which involves zero stretching energy. If denotes the 

plane-stress metric then 

E=l£ - U)isf6 - S6y) dS = Eps, 

whereas if h^^ is the second quadratic form that min- 
imizes the Willmore functional (subject to the satis- 
faction of the Gauss-Mainardi-Codazzi equations with 
gap = gap) then the energy reduces into 

kJ s 

Clearly, if the Willmore energy is lower than the plane- 
stress energy then the plane-stress solution is imstable. 
This provides a lower bound for the buckling threshold, 
known as the crossover point. 

Below this thickness, we expect the solution to approach 
an isometric immersion of minimum energy, thus the 
energy should be close to Ewf- Obviously, in order to 
evaluate the crossover point one needs to know the min- 
imizer of the Willmore functional, which may be highly 
non-trivial (it requires in particular the solution of an 
isometric immersion problem). 

Examples Consider the elliptic and flat geometries 
(11) and (12), for which the isometric immersion that 
minimizes the Willmore functional is explicitly known. 
It is a stirface of revolution, 

fir, 6) = {(Pir) cos 6, ct>{r) sin 6, t/^(r)). (18) 

The corresponding metric is 



hence the isometric immersion satisfies 

(f) = ^ and {(p'f + {ff = 1. 

(For the hyperbolic metric (13) 0'(r) > 1 hence there is 

no axisymmetric isometric immersion.) 

For the elliptic geometry (11) with the same parame- 
ters as above we find 

Eps = 0.0163 and Ewf = 0.3609, 

from which we get t^-.o. = 0.2125, which is lower than 
tb = 0.367 by about 40%. In contrast, we obtain for the 
flat metric (12) 

Eps = 0.0580 and Ewf = 75.64, 

from which we get f^.o. = 0.0277, which is lower than 
tb = 0.387 by more than an order of magnitude. This 
demonstrates that in certain cases the crossover point 
may provide a very poor estimate of the buckling thresh- 
old. The reason why the discrepancy between tc.o. and 
tb may be large is that the buckling transition is a prop- 
erty intrinsic to the plane-stress solution, not to isometric 
immersions. 



D. Bifurcation analysis 

In this section we analyze the nattire of the buck- 
ling transition. At t = tj, the plane-stress solution is 
marginally stable. In particular, there exists a non-trivial 
perturbation that does not change the elastic energy up 
to terms that are quadratic in v, w. Specifically, 

5Ef°\v) = if and only if i; = 0, 
and there exists a it) such that 

6Ef\w) + tl6Ef\w) 

changes sign att = t\,. In fact, as 6Eg can be identified 
as an inner-product (cf. (16)), it follows that for every w, 

{w,'Kw) = Q. (19) 

Since w is determined up to both additive and multi- 
plicative constants, we will define w to have zero mean 
and be normalized, \\'w\\2 = 1, where || ■ II2 is the norm. 

For plate thickness below t}, the flat configuration is 
linearly unstable. (Note however that the plane-stress 
solution is a critical point of the energy functional for all 
values of t; it only ceases to be a local minimum at t},). 
The loss of stability of the flat solution is due to a biftir- 
cation. A branch of stationary solutions with non-zero 
bending content merges with the plane-stress solution at 
t = t},. The bifurcation is called super-critical (forward) if 
the branch of buckled solutions exists for t < t;,, in which 
case, as predicted by bifurcation theory [26], the buckled 
solutions near t], are linearly stable. For a super-critical 
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bifurcation the transition from the plane-stress solution 
to the buckled solution, as t decreases below tt, is contin- 
uous. The bifurcation is called sub-critical (backward) 
if the branch of buckled solutions exists for t > t},. In 
this case, the buckled solutions near tj, are unstable (this 
branch of solutions becomes stable after it turns back). A 
transition to linearly stable solutions occiirs discontinu- 
ously att = t},. In particular, discontinuous bifurcations 
exhibit hysteresis. 

To analyze the bifurcation we need to study the be- 
havior of the energy functional in the vicinity of the bi- 
furcation threshold. Since the plane-stress solution is 
marginally stable at tt, terms that are of higher order in 
V, w must be taken into account. 

Let gap and Sap be the plane-stress metric and stress, 
and set - - e with e > a small parameter. That 
is, we consider plate thicknesses just below the buck- 
ling threshold. By the above discussion, the bifurcation 
is super-critical if for small e the energy fimctional has 
a local minimum for a non-trivial perturbation whose 
magnitude vanishes as e i 0. If the bifurcation is sub- 
critical, then the stable solution for e > does not con- 
verge to the plane-stress solution as e J, 0. Our working 
h5^othesis is that the bifurcation is super-critical. The 
analysis will prove us wrong if this is not the case. 

Set once again 5f = dyf + wJ^. Substituting the 
variations (B6), (B16) in stretching and bending content 
densities, the variation in total energy takes the form 



(5£ 



+ {tl - e) {bEf^\w) + bE^l'^\v,w) + bEf^\w)] 



where 5Ef°\v), 6E'g 
and 



and 6Ef^\w) are given by (14) 



We are seeking the perturbation that minimizes the en- 
ergy variation for small e > 0. Since we expect, to lead- 
ing order, the minimizer to be proportional to w (the 
least stable out-of-plane mode at t},) with a pre-factor 
that vanishes as e — > 0, we expand the minimizer in a 
power series in e, whose first terms are 



i;>' = cVeP-l-... 
w = cwe^ + We'' + 



(20) 



where the exponents p, q, r and the constant c are yet to be 
determined. Substituting this expansion into the energy 
variation we get 

[6Ef"\v)] e^^ 
c" [dEf'-\v,w) + tl dE^l''-\v,w)\ e^^^ 
[6Ef\w)^tl6Ef'\w)\ 



+ c 



+ c 



+ 0{e^', e^P, e"^, e^P+'J, eP+S"?, e^p+i^ ^4,+!^ g.p+2<?+i) 

The first term on the right hand side is the quadratic 
out-of-plane term, which is, as expected, negative for 
e > 0. For a super-critical bifurcation. It is balanced by 
the quartic term, from which we infer that lq+1 = Aq, i.e., 
q = 1/2. Since the term that is quadratic inu is positive, 
it follows that 2p > 2q + 1 = 2. We may then set p = 1, 
with the possible outcome that we obtain ^ = (i.e., that 
the V terms are sub-dominant). 

We proceed to minimize this expression (with all four 
terms of order e^) with respect to the in-plane perturba- 
tion V and the constant c. Note that imder the normaliza- 
tion choice in (20) the mirdinizing v does not depend on 
c, since the two i;-dependent terms are proportional to 
c^. The 0-dependent terms consist of a positive-definite 
quadratic term and a linear term, which guarantees the 
existence of a non-trivial minimizer (and in particular 
confirms that p = 1). Once v has been determined, a 
minimizing c exists if and only if the sum of the terms 
proportional to are positive. Then, 



dE^iw) 



6Ef^\v) + [6Ef'\v, w) + tl 6£<i'2>(f), w)] + [dEf'\w) + tl bEf'\w)\ ' 



Recall that ec is, to leading order in e, the L^-norm of the 
out-of-plane deflection w. If the denominator is nega- 
tive, then the bifurcation is sub-critical and the branch 
of stable buckled solutions cannot be foimd by a local 
analysis about the plane-stress solution. 

We calculated c for both the elliptic and flat geometries; 



recall that in both cases the critical mode is axisymmetric, 
n = 0. For the elliptic geometry the bifurcation was 
found to be super-critical for the whole possible range 
of curvatiires K (for large enough K the surface is no 
longer an embedding, as the sphere closes upon itself). 
For the flat geometry a transition from super-critical to 
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FIG. 7: The coefficient c versus the Gaussian curvature K of the 
elliptic geometry (a) and the parameter a of the flat geometry 
(b). 

sub-critical bifurcations was found: the bifurcation is 
super-critical for a in the range 0.58 < a < 1, and sub- 
critical for a < 0.58. 

In Figure 7 we show the value of c versus the Gaus- 
sian curvature K of the elliptic geometry (left) and the 
parameter a of the flat geometry (right). Note that at the 
transition point from super-critical to sub-critical bifur- 
cation, a ~ 0.585, the coefficient c diverges. 

V. BOUNDARY LAYERS IN VERY THIN PLATES 

It was shown in Section III that provided that a limit 
equilibrium configuration as f — > exists, it is given 
by the isometric immersion that minimizes the Will- 
more functional. How a sequence of equilibrium con- 
figurations approaches the Willmore isometry is non- 
trivial. The convergence is in the Sobolev space ^ — 
the space of surfaces with square integrable second 
(weak) derivatives [27]. This guarantees (by the Sobolev 
embedding theorem) uniform convergence in the space 
of once-differentiable embeddings, but not in the space 
of twice-differentiable embeddings. In other words, sec- 
ond derivatives may not converge uniformly. 

Almost a hundred years ago (see [28] and references 
therein), it was observed that thin elastic bodies may ex- 
hibit boundary layers, which interpolate between a state 
of minimum stretching content in the bulk and the zero 
normal traction and zero bending moment conditions at 
the botmdary. Such boundary layers also occur in non- 
Euclidean plates, and turn out to dominate the deviation 
from an isometry as t — > 0. 

Generally speaking, a large thickness implies a bend- 
ing energy-dominated configuration (i.e., close to flat), 
whereas a small thickness implies a stretching energy- 
dominated configuration (i.e., close to an isometry). 
Whether a thickness t is to be considered as "large" or 
"small" is determined by comparison with the shortest 
lengthscale of the problem, which may vary with posi- 
tion. For every finite t there exists a distance from the 
boundary, {, with respect to which t cannot be considered 
small. As a result, we expect bending energy-dominated 
behavior in a strip of thickness { near the boundary. 

We start with a scaling argument. Let ha^ be the second 



fundamental form of an isometric immersion f{x^,x^) 
that minimizes the Willmore functional (the metric is of 
course equal to the reference metric ga^ = ga^). From 
the point of view of the bending energy it would be 
favorable to have a flat surface, ha^ = 0, however the 
second fundamental form cannot be modified without a 
modification of the metric, as the two must satisfy the 
Gauss-Mainardi-Codazzi equations. In particular, the 
Gaussian curvature is an isometric invariant. From the 
analysis in Appendix B, assuming that the perturbations do 
not involve small-scale features, we see that the variation in 
stretching content is quadratic in the perturbation fields 
V, w, whereas the variation in bending content is linear 
in V, w, i.e., 

5E~ 0{v^,w^) + t^ 0{v,w). 

Since equilibrium is obtained by a balance of the neg- 
ative bending contribution and the positive stretching 
contribution, then v,w ~ 0{t^), and dE ~ 0{t'^). 

These are however bulk considerations, where energy 
balance is considered uniformly over the surface. The 
question is whether the total elastic energy can be re- 
duced by a larger, yet local change in the bending content 
density. This is not possible inside the domain, because 
even a local change in the second fundamental form in- 
volves a non-local change in the metric, hence the gain 
in stretching energy exceeds the loss in bending energy. 
The situation may however be different at the boundary. 

In a bending-dominated region we expect the curva- 
tures to deviate from the curvatures associated with the 
Willmore isometry by 0(1). As curvatures relate to the 
metric through two differentiations, such a deviation 
over a strip of width tP, induces a metric deviation of 
0{t^). Thus, the variation in total energy is of order 

6E ~ O(f^P) + 1^ 0(1), 

which yields, p = 1/2. 

To make this into a rigorous argument, suppose that 
v = 0{ti) and w = 0(f), where q,r > 0, both varying over 
a boundary layer whose width scales like f, where p >0. 
Inside the boundary layer, the variation in stretching 
energy density (B6) is dominated by terms of order 

6Ws = 0(t2'!-2P,p,f'?+'-P,f?+2'-3p^t3'-2p^f4r-4p^^ 

Note that this contiibution is always positive. On the 
other hand, the variation in bending energy density, 
(B16), which can become negative, is dominated for 
small t by terms of order 

^^6l^;B = 0(^'^-P+^r2p+2). 

The exponents q, r, p are determined such to maximize 
the change in (negative) bending energy, without the 
(positive) stretching energy becoming dominant. That 
is, if we define 

es = min(2(j -2p,2r,q + r -p,q + 2r- 3p, 3r - 2p, 4r - 4p) 
Cb = mm{q -p + 2,r -2p + 2), 
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then we need to choose q, r, p such to minimize eg subject 
to the constraint that eg > £b- It can be shown that the 
optimal choice satisfies r = 1, p = 1/2, and q > 3/2. 
That is, the width of the boundary layer is expected 
to scale like the square root of the plate thickness. To 
minimize the gain in stretching content, we expect a 
balance between the V„c'^ and w terms, which yields 
q = 3/2. 

Assuming these scaling exponents, we study the struc- 
ture of the boundary layer. We consider, as before, a 
perturbation, which we decompose as 

Consider now a local parametrization / : [0, X 
[0,{^) —> M.^ of the surface, such that the parametric 
line = coincides with a boundary of the surface, 
with the positive axis inside the sample. Moreover, 
the parametrization of the unperturbed surface is semi- 
geodesic, i.e., gn = 1, gi2 = gii - and g22 = f/)^; such 
a parametrization is always possible. One may also set 
g22 = 1 along the boundary (see Figure 8). 

Since we expect a boundary layer of size y/i, we stretch 
the positive axis accordingly by introducing a rescaled 
coordinate, 

Vf 

and rescale the perturbations v'',zv, such that the new 
variables and their derivatives are all of order one, 

W{X\x^) = (VfXi,x2) 

\' (21) 

By setting, say, i'^ ~ 0(f^^*) we have a situation where, as 
t — > 0, the local coordinates {x^,x^) parametrize a shrink- 
ing annulus which converges to the boundary, whereas 
in the rescaled coordinates, (X^,x^), the range of X^ in 
the positive direction tends to infinity. We are going to 
show the existence of a perturbation of such structure 
that reduces the total elastic energy. 

We then evaluate the variation in energy content den- 
sities inside the boundary layer, i.e., at points = Vfx^, 
with X^ ~ 0(1)- To leading order, the unperturbed met- 
ric and the Christoffel symbols are given by their values 
at the boundary, and covariant derivatives coincide with 
partial derivatives. Since gap = gafi we also have 

A"''^" = Y^5ati5y6 + \ [bayb^b + dasbpy) + 0(fi/2). (22) 

Substituting the rescaled variables (21) into the vari- 
ations (B6) and (B16) in stretching and bending content 




FIG. 8: (Color online) Local parametrization of an annulus 
bounded by a boundary of the domain. 

densities, we get 

6ws/t^ = ^A'^'\diVt'){diV'') + ^A^PrXphy.W^ 

- A^p^%6{diVi^)w + iyi""(,9iyi)(^i w)2 

- lA"l^^^ha.{diWfW + iyi""(<9i W)* + 0{t^'% 
2 8 

and 

6wB = ^A"^''h„ii{d,d,W) + ^yi""(<9i5iW)2 + 0{t^'^). 

Substituting also expression (22) for the elastic tensor, 
we end up with 

6ws/t^ = \ [(<9iy2) - 2/zi2W]' + l/z2^W2 

+ g^Y^:^ [2(^1 1^') - 2(^11 + ?!22)w + (5i w)2]' 

+ i \l{dxV^) - 2hnW + {diWff + 0{t^'^) 

bzvB = ^^(^^ [2 (hn + vh22) (^i^i W) + {did.Wf] 
+ 0(fi/2). 

The variation in stretching content density is a sum of 
squares. The optimal choice of is the one that makes 
the first term vanish, namely, 

(diV^) = 2Whi2. 

Similarly, the optimal choice for satisfies 

2{diV^) + {diWf = 2 (/ill + vh22) H 
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so that the variational problem reduces into one for W 
only. 



y',y2 optimal 



Omitting the 0{t^^^) terms, the resulting Euler-Lagrange 
equations are 



5(Xi)4 

with boundary conditions 



+ 12(l-v2)/i^2W = 0, 



= - (/ill + v/122) , and 



= 0, 



,9(Xi)2 ' 5(Xi)3 
at = 0, and W(oo, x^) = 0. The solution is 

WiX\2) = _^!ll±^e-AX^ (cosAXi - sinAXi), (23) 

where 

A = [3(l-v^)/iy' 



ll/4 



Reverting to the original unsealed units, {x^, x^), the out- 
of-plane perturbation w{x^,x^) exhibits a boundary layer 
whose width ^b.l. is 



Substituting the as5anptotic boundary layer profile 
(23) in the leading order expressions for bwsff' and bwB, 
we get 



, = %±^e-^^^^ (cos AXi - sin AX^f 



24(1 - v) 

(/Zll + r/;22)^ _2AXi 

Swr = —e 



24(1 - v) 
X (cos AX^ + sin AX^ - 2 e^^') . 



(cosAX^ +sinAXi) 



(24) 



At the boundary, to leading order, curvature in the 
direction normal to the boundary is given by 

/ill + d\diw = -vhzz. 

This leads to the vanishing of the bending moment at 
the boundary, i.e., = 0, which is one of the bound- 
ary conditions associated with the energy minimization 
variational problem [1]. Note that the normal bending 
moment along the bovindary is proportional to /in + v/122. 
If V = and the bending minimizing isometry satisfies 
hu = 0, then no correction is needed in order to satisfy 
the boundary conditions, and we expect no boundary 
layer to develop. This fact is manifested by the van- 
ishing of W in (23). For v + Q, however, the curvatures 
normal and tangent to the boundary have opposite signs. 



which means that the surface is h3^erbolic in the vicin- 
ity of the boundary. Finally, our analysis breaks down in 
the event that /122 = 0, in which case a boundary layer of 
different nature may emerge. 

To validate our results, we plot in Figure 9 the rescaled 
deviations in stretching content density, bwslf", and the 
deviations in bending content density, bw^, versus the 

rescaled coordinate {Vmax - Vf, for the elliptic geom- 
etry (11), with the same parameters as in previous sec- 
tions, namely, rmin = 0.1, r^ax = 1-1/ i*^ = 1/ and v = 0. 
The results were obtained by minimizing the full energy 
functional. We display results for t = 0.01, 0.005, 0.002, 
and 0.001. As expected, the four rescaled curves approx- 
imately coincide. These niraiericaUy computed curves 
are compared with our asymptotic expressions (24). 



VI. DISCUSSION 

A theory of non-Euclidean plates, applicable to thin 
elastic sheets that do not have a stress-free rest config- 
uration, has recently been proposed in [1]. This paper 
provides a first mathematical analysis of this model. Two 
different limits of such plates are analyzed: (i) the buck- 
ling transition, and (ii) the occurrence of boundary layers 
in the limit where the plate thickness tends to zero, and 
the configuration is expected to converge to the isometric 
immersion that minimizes the Willmore functional. 

We proved a general result, whereby any non- 
Euclidean plate that does not have a flat stress-free 
configuration (i.e., whose reference metric has non-zero 
Gaussian curvature), buckles if the plate is sufficiently 
thin. The transition from flat to buckled equilibria may 
be either continuous, or discontinuous, depending on 
the particular system. Instances of both tj^es have been 
observed. 

We showed that in the thin plate regime, the dominant 
deviation from the Willmore isometry is governed by 
a bending-dominated boundary layer, whose structure 
was calculated using a boundary layer analysis. In par- 
ticular, the width of this boundary layer is determined 
by both the plate thickness and the tangential curvature 
at the boundary of the Willmore isometry — it scales like 
the square root of their product. It would be of interest 
to observe the occurrence of such boundary layers in ex- 
periments, e.g, in the thermo-responsive gels studied in 
[4], in order to further validate the model as describing 
thin elastic sheets with no stress-free configuration. 

The reported results are of high relevance to shape for- 
mation in growing biological tissues . Spontaneous buck- 
ling and wrinkling of non-Euclidean sheets was sug- 
gested as a mechanism for shaping leaves that are free of 
external confinement [2, 3, 23, 29]. These studies were ei- 
ther qualitative, or assumed the limit of an infinitely-thin 
sheet. As such, their results are relevant only to selected 
cases. Recent studies, suggest and demonstrate that the 
mechanical stress field can lead to differentiation of cells 
[30-32], and might act as a regulator of tissue growth [33] . 
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FIG. 9: (Color online) Structure of the boundary layer near 
the outer boundary, r = 1.1, for the elliptic geometry (11) and 
the same parameters as in previous sections. The top figure 
(a) shows the rescaled deviation in stretching content density, 

bwsif versus the rescaled coordinate (r^ax - '')/ Vt, calculated 
for four different thicknesses, h = 0.01, 0.005, 0.002, 0.001 (sym- 
bols). The solid (blue) line is the asymptotic result (24). The 
inset shows the unsealed stretching content density deviations 
versus the unsealed coordinate - r. Note the logarithmic 
y-scale. The bottom figure (b) shows the bending content den- 
sity, bwB, versus the rescaled coordinate for the same values of 
the thickness. The inset shows the unsealed results. 



These studies emphasize the need to know the mechan- 
ical state of a tissue of finite thickness, which undergoes 
differential growth. Our calculations of buckling thresh- 
old, as well as the pre- and post-buckling stress distribu- 
tion within plates, can be integrated into the biological 
picture as inputs that can affect its future evolution. In 
particular, the existence of a boundary layer and its scal- 
ing predict thickness-dependent, localized effects near 
the boundaries of unconstrained growing tissues. 
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APPENDIX A: THE t LIMIT 

The t ^ limit can be approached in two way. The 
first possibility is to depart from the three-dimensional 
model, i.e., the energy functional (1), and analyze the 
limit of the energy minimizers (or approximate energy 
minimizers) as f ^ 0. To this end, one would hope 
to be able to use the analytical techniques based on T- 
convergence developed in [34-36]. There is however 
an obstacle that prevents the direct application of the 
abovementioned techniques to the present context: we 
do not have a reference configuration with respect to 
which deviations can be measured. Indeed, the analysis 
in [35] relies heavily on a so-called rigidity property that 
estimates the distance of the displacement from a rigid 
rotation in terms of an integral over local distances from 
rotations. 

The second alternative is to depart from the two- 
dimensional model, i.e., the energy functional (6). For 
reasons to be clarified below, we work with an energy 
functional rescaled with the thickness square. 



E _ 1 

f2 " f2 



(Al) 



where the notation Ft makes the f-dependence explicit. 
Clearly, for every fixed t the functionals £ and Ff have the 
same minimizers. We view Ft as a one-parameter family 
of functionals, defined on the Sobolev space W^'^(S;]R^), 
of surfaces whose second (weak) derivatives are in L^(§). 
Since we view every two configurations that differ by a 
rigid motion as identical, the space of immersions is in 
fact the quotient space W^'^(§; ]R^) modulo rigid motions. 
We denote by Ff[/], £s[/]/ and £b[/] the functionals Ff, 
£s, and Eb evaluated at a configuration / = f{x^,x^). We 
will also denote by 

et = mf[FAf] ■■ /€ W2'2(§;]R3)} 

the f-dependent greatest lower bounds on the energy. 

The two-dimensional elastic problem formulated in 
the Section II assumes the existence of a family of mini- 
mizers, /j*, such that Ff[/f*] - et- Even if such minimizers 
do not exist, we can always construct a family of approx- 
imate minimizers, /*, satisfying, 

lim(Ft[/t*]-et) = 0. 
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Suppose now that the two-dimensional reference met- 
ric gap assumes an isometric immersion f = f with finite 
bending content. Then for every t, 

et < F,[f] = EBlf], 

i.e., the greatest lower boimds on the energy are imi- 
formly boimded. In particular, it follows that 

limEs[/;]=0, 

which means that the metrics g{ff) associated with the 
family of approximate minunizers converge, in least- 
square norm, to the reference metric g. 

The mean-square convergence of the family of met- 
rics, as t ^ 0, does not guarantee that the family of 
(possibly approximately) minimizing configurations has 
a limit (modulo rigid motions). If, however, such a limit 
does exist, then we show that this limit is an isometric 
immersion that minimizes the bending content, i.e., the 
Willmore fimctional. Specifically, 

Let g be a reference metric that assumes a W^'^ 
isometric immersion. Let be a family of ap- 
proximate minimizers of the functionals Ft- If the 
family fl (weakly) converges in W^'^(§;]R^), as 
t ^ 0, to a limit /*, then f* is a configuration 
that minimizes the Willmore functional over all 
isometric immersions of the reference metric g. 

To prove this theorem we construct a "limit func- 
tional", 

\oo otherwise, 

and show that the fiinctionals Ft F-converge to Fq, as 
t — > 0, with respect to the weak W^'^ topology. It then 
follows that every converging sequence of approximate 
minunizers of F; converges to a minimizer of Fq [37]. 

To show that Ft T-converges to Fq we need to show 
that: 

1. Lower-semicontinuity: for every sequence /( that 
converges to a configuration / (in the weak W^'^ 
topology), 

limmf Fi[/,] > Fo[/]. 

2. Recovery sequence: for every / e W^'^ there exists 
a sequence ft that weakly converges to / for which 

liminf Ff[/f] = Fo[/]. 

To prove the lower-semicontinuity property we note 
that the weak W^'^ convergence of ft to / implies the 

weak convergence of the corresponding metrics, gift) — > 
gif), in the weak W^'^ topology, which by the Sobolev 



embedding theorem [27] implies the convergence of the 
metrics in the strong C° topology, i.e., uniform conver- 
gence. It follows at once that the corresponding family 
of second fimdamental forms weakly converges in to 
the second fundamental form of /, h{ft) — > h{f). Since 
the bending content is equivalent to an norm of the 
second fimdamental form, it follows at once that 

liminf EB[/i] > Eb[/]. 

t->o 

Now either g{f) = g, in which case 

limF,[/(] > liminf Eb[/,] > Eb[/] = Fo[/], 

or g{f) g,in which case 

liminfFt[/f] = oo = Fo[/]. 

To prove the existence of a recovery sequence we take, 
given /, the constant sequence, ft = f. If g{f) = g then 

limFi[/f]=EB[/]=Fo[/]. 

If, however, g(f) + g then 

limFi[/f] = cx, = Fo[/]. 

f->0 

Comments 

1. The assumption that the reference metric can be 
embedded isometricaUy with finite bending is by 
no means trivial. The Nash-Kuiper embedding 
theorem only guarantees the existence of a em- 
bedding. Embeddings of class W^'^ have been 
shown to exist under additional assirmptions (see, 
e.g., [38]), however there is no general existence 
proof for arbitrary metrics. 

2. We use the weak W^'^ topology because we aim to 
eventually prove that every family of approximate 
minimizers has a converging subsequence (imply- 
ing that the limit is a minimizer of the Willmore 
functional). Such a compactness result cannot pos- 
sible hold in the strong W^'^ topology. 

3. A side-result of the above theorem is that the Will- 
more functional has a minimizer (although not nec- 
essarily unique), even if the functionals Ft do not 
have minimizers. 



APPENDIX B: PERTURBATION ANALYSIS 

Consider a sufficiently regular surface f{y} ,■)?■). Any 
small enough perturbation can be decomposed into a 
sum of in-plane and out-of-plane displacements, 

bf = vy dyf + w'H, (Bl) 
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where N is the unit vector normal to the surface. Given 
such a perturbation we are going to calculate the varia- 
tion in the elastic energy. 

In order to retain the tensorial nature of the problem 
(coordinate invariance), we need to only utilize covari- 
ant differentiation. To do so we need to specify a metric 
with respect to which the Christoffel symbols are de- 
fined. It turns out that choosing the (natural) induced 
metric on the surface, yields the most compact form for 
the variation in energy. 

We recall the definitions of the covarian.t derivatives 
of a scalar field W, a contravariant vector W, a covariant 
vectors Vy, and a mixed tensor Ty, 

V„W = daW 

« (B2) 

'7aVy=daVy-T^^yVp 
Vaiy - Cai^ -t- l^^i^ layl^. 

Note that and Vjs commute only when operating on 
scalars. As operators on higher rank tensors their com- 
mutator is nonzero, and relates to the Gaussian ciirvature 
of the surface. 

To calculate the variation in energy we need to calcu- 
late the variation in the two fundamental forms. For this 
we use the Gauss- Weingarten equations, 

dad^f = Tl^dyf + hapS and daH = -T^ d^f, 

where = {g~^)^^hya. (Note that by our definitions of 
index raising, = h^. only if ga^ = gap.) It follows that 
for a vector in in the form 

its derivative is given by 

dpv = (V^a« - bTfj daf + (Vpb + a^hap) (B3) 

1. Variation in stretching content 

Differentiating (Bl) and using (B3) we get 

da^f = [i^^aVn - ^ Tl] dyf + [{V,W) + vy hay] K (B4) 

from which we derive the variation in the metiic, 

6gaj3 = daf ■ dfibf + dabf ■ dpf + da6f ■ dpdf 
= gtSyiyaV^) + gay{ypV^) - 2whap 

+ [{VaU>) + vyhay] [{V,,W) + V%6] ^^^^ 
+ [{Vavy) - wTll gy, [(Vpv') - WT^] . 



Substituting into (5) we obtain the variation in stretching 
content density, 

6ws = \s"^dgap + l-A'^^^'^gap^gyi 

= 6w'-^^\v) + 6wf^\w) + 6wf'°\v) + 6wf^\w) 
+ 6wf^\v,w) + 6wf^\v,w) + 6wf'^\w) + 6wf^\w) 

(B6) 

where the various 6wg''^ represent terms of different or- 
ders in V and w, 

6wf^\zv) = -s"l^haiizv 

dwfyv) = ^s"fihayhpr,vyv^ + ^s"^g^e(V„i;>')(V^i;^) 
+ lA"^'"gp,g5eiyaV'^)iyyV') 

+ ^A'^i'y^phyaw^ 

6wf^\v,w) = s"%yVy{VaW) - s"^ha^iWpV'')w 
- A"^y%,gpr,iWaV^)w 

6wf^\vM = ^A^^ygp^iVaV^) [{yywmw) + Hyew^] 
-A"^y%aph^,v^w{Vyw) 

+ A"^y%aphy{yyv^') 
dwf^\w) = -^A"i^y%p [{yyw){y5w) + Hy^w^] w 

bwY\w) = ^A^l'y' [(V„it;)(V^u;) + Hapw^] 

X [(VyM')(V6ll') + Hy5W^] , 

where we have used the syrranetiy of s"^ and A^^y^, and 
intioduced the following new S5anmetric tensor. 

Hap = Tlhyp = {g-')y%ahyp, 

which is known as the third quadratic form. 

Perturbation of a flat surface When the unperturbed 
surface is flat then hap = Hap = 0, which simplifies (B6) 
considerably. In particular, all the terms that are odd 
functions of the out-of-plane perturbation w vanish. 

Perturbation of an isometric immersion If the im- 

perturbed surface is an isometric immersion, gap = gap, 
then s"^ = 0, which implies that the lowest-order non- 
vanishing terms in (B6) are 

6WS = ^A"Py' {gp^iVaV'^) - hapw) (gyei^ 6V') - hy,w) 
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Note that it is explicitly assumed here that derivatives of 
the deviations v, w are of the same order of magnitude as 
the deviation themselves. This assumption breaks down 
in the presence of small scale features, such as boimdary 
layers. 



2. Variation in bending content 

To calculate the variation in the second quadratic form 
we start with 

hap = dadfif ■ N = -dpf ■ da'N. 

from which follows that 

6Kp = -dp6f ■ d^J^ - dpf ■ (9„6N - diidf ■ dadn. (B7) 

The first term follows directly from (B4) and the Wein- 
garten equation, 

-dp5f ■ 5„N = hayiV^v^) - wHa^. (B8) 

To calculate the second and third terms, we need to 
express the pertiirbation 6N in the unit normal vector. 
For that, we use the facts that 6(N ■ N) = 6{daf ■ N) = 0, 
hence, 

2N • 6N + 6N • 6N = dadf • N + 5«/ • + da5f • 6N = 0. 
Setting 

d3s[ = {g-y%dyf + bK (B9) 

this yields a closed set of equations for the three coeffi- 
cients Uy, b, 



fl„ = -(1 + b) [{y„w) + vyKy] - [{Vy) - wT^] 
h = -\{g-y%ap-\b\ 
Applying (B3) on (B9) we get 



(BIO) 



<?a6N = [(g-^y'Vad, - bTl] dyf + [V„b + Tiae] K 

(Bll) 

where we have used the fact that the covariant derivative 
of gap and its inverse vanishes. It follows at once that the 
second term in (B7) is given by 



-dpf ■ da6J< = -{Vattp) + bhaf}. 



(B12) 



The third term in (B7) is obtained by combining (Bll) 
and (B4), 



-dpbf ■ dab-N = - [{Vpvy) - wTj^] [{Vaay) - bhay] 



[(VpW) + vyhpy][(\/ab) + '!%]. 



(B13) 



In remains to combine (B8), (B12) and (B13) to get 

- (Vafljs) + bhafl 

- [(Vpvy) - WTI\ [(Vaay) - bhay] ^^^^^ 

-[(Vpit;) + i;''V][(V«b) + T%] 

So far all the relations are exact, i.e., no assumptions 
have been made about the smallness of the perturbation, 
other than the ability to decompose it in the form (Bl). 
Equations (BIO) are a set of three quadratic equations for 
fly, b, which we may solve by successive approximations. 

Perturbation of a flat surface When the siirface is flat, 
hafi = 0, (BIO) and (B14) reduce into 

h = -\ig-y^ayap-\b\ 

and 

bhap = -(V„flp) - iVpvy)iVaay) - iVpw)iVab) 
Solving for Ua, b by successive approximation, we get 

aa = -(V«a;) + {Vpw){Vav'') 

+ lig~yHyaW)iypW)iVyW) + Oiv\ VW\W^) 

Putting it all together we obtain a simple expression for 
the variation of the second form, 

5hap = (V^V^jo;) - (Vj,a;)(V„V^i;>') 

- ^(^"')'"(VaV^a;)(VyM;)(V6w) + 0(v\vw\w^). 

(B15) 

We then substitute (B15) into (5) to calculate the variation 
in bending content density, 

(B16) 

where 

6wf^\w) = ^A"Py\WpWaZV){V,VyW) 

6wf^\v,w) = -^A''i'y'{y^Vavy){yyw){VyV6iv) 

6wf'\w) = -^A"^y'ig-'riWr,w)iWeW)iWpWaW)iW5'^yW). 
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